%--------------------------------------------------------------------------
% CDS PRICING
%--------------------------------------------------------------------------

% =========================================================================
% default probabilities under Q
% =========================================================================

[Cheb_nodes,~] = qnwcheb(nE,-1,1);
k_vec       = repmat(k',[nS,1]);
k_vec       = k_vec(:); 
kb_vec      = repmat(k_bar,[nK,1]);
kb_vec      = kb_vec(:); % [nS*nK,1]

mu_vec      = repmat(mu_XQ,[nK,1]);
mu_vec      = mu_vec(:); % [nS*nK,1]
sig_vec     = repmat(sig_X,[nK,1]);
sig_vec     = sig_vec(:); % [nS*nK,1]
Q_mat       = repmat(Q,[nK,1]); % [nS*nK,nS]

default     = reshape( k' <= k_def ,[nK*nS,1]);
issuance    = reshape( k' >= k_iss ,[nK*nS,1]);

Amat        = (log(k_def)' - mu_vec - log(k_vec) )./sig_vec;
Amat_iss    = (log(k_def)' - mu_vec - log(kb_vec))./sig_vec;
Amat(issuance,:) = Amat_iss(issuance,:);

q1          = reshape(~default.*sum( normcdf(Amat) .* Q_mat ,2),[nS,nK]); % [nS,nK]
q           = zeros(nS,nK,Tmax_CDS);
q(:,:,1)    = q1;

for t=2:Tmax_CDS % maturity
    for is=1:nS % future state
        a       = max(-5,(log(k_def(is)./k_vec) - mu_vec)./sig_vec); 
        b       = 5;     
        Eshocks = (b-a)*(Cheb_nodes'+1)/2 + a;
        temp    = pi*(b-a)/(2*nE);
        weights = temp.*sqrt(1-Cheb_nodes'.^2);
        probE   = normpdf(Eshocks).*weights;
        
        % kappa transition
        kp             = k_vec  .* exp(mu_vec + sig_vec.*Eshocks); %
        kp_issuance    = kb_vec .* exp(mu_vec + sig_vec.*Eshocks); %
        kp(issuance,:) = kp_issuance(issuance,:);
        kp(default,:)  = k_min; % any default prob evaluated at k_min equals zero
        
        q_Fit    = griddedInterpolant(k,q(is,:,t-1)','linear');
        q(:,:,t) = q(:,:,t) + reshape(Q_mat(:,is) .* sum(q_Fit(kp).*probE,2),[nS,nK]); % [nS,nK] %%%
    end
end


% =========================================================================
% loss given default under Q
% =========================================================================


lamD_fit        = griddedInterpolant({1:nS,k},lamD,'linear');

mu_vec          = repmat(mu_XQ,[nK*nS,1]);
sig_vec         = repmat(sig_X,[nK*nS,1]);

k_vec           = repmat(k',[nS,1,nS]);
k_vec           = k_vec(:);

VEC_k_bar       = reshape(repmat(k_bar',[nK*nS,1]),[nK*nS^2,1]);
kb_vec          = repmat(k_bar',[nK*nS,1]);
kb_vec          = kb_vec(:);

VEC_lamD        = reshape(repmat(lamD_fit((1:nS)',k_bar)',[nK*nS,1]),[nK*nS^2,1]);

MAT_Q           = reshape(permute(repmat(Q,[nK,1,nS]),[1,3,2]),[nK*nS^2,nS]);
default         = reshape(repmat(k'<=k_def,[1,1,nS]),[nK*nS^2,1]);
issuance        = reshape(repmat(k'>=k_iss,[1,1,nS]),[nK*nS^2,1]);

Amat            = (log(k_def)' - mu_vec - log(k_vec) )./sig_vec;
Amat_iss        = (log(k_def)' - mu_vec - log(kb_vec))./sig_vec;
Amat(issuance,:) = Amat_iss(issuance,:);

MAT_cdf         = normcdf( Amat-sig_vec );
MAT_Rf          = permute(repmat(Rf(:,1:Tmax_CDS),[1,1,nK,nS]),[1,3,4,2]);
L1              = repmat(q1(:),[nS,1]) - MAT_Q .* MAT_cdf .* k_vec .* exp(mu_vec+sig_vec.^2/2) ./ (VEC_lamD.*VEC_k_bar) * ((1-omega).*lamA);
L1              = max(L1,0);
L1              = reshape(L1,[nS,nK,nS]);

% transition of most recent issuance state (note: predetermined)
markov_state     = repmat((1:nS)',[nK*nS,1]);
xi_s_p           = reshape(repmat(1:nS,[nK*nS,1]),[nK*nS^2,1]);
xi_s_p(issuance) = markov_state(issuance);
xi_s_p_dummy     = zeros(nK*nS^2,nS);
xi_s_p_dummy((xi_s_p-1)*nK*nS^2 + (1:nK*nS^2)') = 1;

% multi-period L
L_q             = zeros(nS,nK,nS,Tmax_CDS);
L_q(:,:,:,1)    = L1;

for t=2:Tmax_CDS % maturity
    
    L_temp = zeros(nS^2*nK,nS); % 2nd dim: future last issuace state
    
    for is=1:nS % future markov state
        
        a       = max(-5,(log(k_def(is)./k_vec) - mu_vec)./sig_vec); 
        b       = 5;        
        Eshocks = (b-a)*(Cheb_nodes'+1)/2 + a;
        temp    = pi*(b-a)/(2*nE);
        weights = temp.*sqrt(1-Cheb_nodes'.^2);
        probE   = normpdf(Eshocks).*weights;
        
        % kappa transition
        kp          = k_vec  .* exp(mu_vec + sig_vec.*Eshocks); %
        kp_issuance = kb_vec .* exp(mu_vec + sig_vec.*Eshocks); %
        kp(issuance,:) = kp_issuance(issuance,:);
        kp(default,:)  = k_min; % any default prob evaluated at k_min equals zero  
        
        for is_s=1:nS % last issuance state state
            L_Fit          = griddedInterpolant(k,squeeze(L_q(is,:,is_s,t-1)),'linear');
            L_temp(:,is_s) = L_temp(:,is_s) + MAT_Q(:,is) .* nansum(L_Fit(kp).*probE,2);
        end
    end
    L_q(:,:,:,t) = reshape(nansum(L_temp.*xi_s_p_dummy,2),[nS,nK,nS]);
end

% CDS rate
cf_PS   = cumsum(L_q./MAT_Rf,4);
cf_PB   = cumsum(permute(repmat(1-cumsum(q,3),[1,1,1,nS]),[1,2,4,3])./MAT_Rf,4);
CDS     = 12 * cf_PS./cf_PB;
CDS(CDS==0) = NaN;

% Loss given default
LGD_q   = L_q ./ permute(repmat(q,[1,1,1,nS]),[1,2,4,3]);
PD_q    = cumsum(q,3);
PD_q(PD_q==0) = NaN;